Evaluation method of surrounding rock stability: Failure approach index theory of strain limit analysis for engineering applications

In general, the ultimate parameter selection method of the failure approach index theory among the three-dimensional problems in geotechnical engineering is unclear in theory, and the symbol convention of the failure approach index in engineering calculation is contrary to the stipulation of the numerical simulation software. Hence, the values of the ultimate plastic shear strain are difficult to determine. To solve this problem, the criterion of positive tension and negative compression and the sequence of the principal stress σ1 ≤ σ2 ≤ σ3 are defined in this paper, and the expression of Mohr–Coulomb yield approach index id deduced. Under the condition of the principal strain sequence ε1 ≤ ε2 ≤ ε3, the formula of the ultimate shear strain is derived using the method of the ultimate strain analysis so as to obtain the simple expression and calculation method of the ultimate plastic shear strain, which has provided the calculation parameters for the three-dimensional ultimate plastic shear strain in the Mohr–Coulomb strain softening model and improved the failure approach index theory. In the light of the aforementioned theory, the ultimate strains of cubic concrete specimens are analyzed, and the obtained ultimate strain values are found consistent with previous research findings, which verifies the correctness and reliability of the ultimate strain analysis method. In addition, it is applied to the quantitative elastic–plastic failure analysis of the section coal pillar in Hengjin coal industry for determining its reasonable retainment width. Consequently, the research results can be embraced as the theoretical basis for the stability analysis of geotechnical materials and exhibits engineering application potential.


Introduction
In recent years, projects like transportation, water and oil transportation pipelines, deep underground disposal of nuclear waste, deep coal mining and, etc. are continuously increasing. Thus, the stability of surrounding rock in underground excavation space must be solved urgently. At this point, the surrounding rock stability evaluation has become a heated issue in academic and engineering circles, and the surrounding rock stability analysis of underground engineering is an important method for comprehensive stability evaluation, as well as an important basis for safety design and construction [1][2][3][4]. Furthermore, the failure approach index theory has been adopted by scholars to analyze the stability of surrounding rock. In this way, the concept of the failure approach index proposed by Zhang et al. [5] can effectively evaluate the risk degree of surrounding rock in the elastic stage and the damage degree in the plastic phrase, which is based on the Mohr-Coulomb yield approach index. Its calculation formula can be expressed as: Here, FAI refers to the failure approach index; ω represents the complementary parameter of the yield approach index (YAI); and FD signifies the failure degree and is represented by Here, γ p denotes the plastic shear strain, while g p f represents the ultimate plastic shear strain. When 0 � FAI < 1.0, the surrounding rock is in the elastic stress state; when 1.0 � FAI < 2.0 and FAI � 2.0, it is in the plastic yield and failure state, respectively.
Above all, the failure approach index can accurately and quantitatively evaluate the stability degree of rock mass in each area of surrounding rock and express the location and range of the failure zone, damage zone, and disturbance zone. Moreover, to predict the influence of excavation disturbance on the stability of rock mass in the future, the evolution law of the stability state of rock mass in each zone along with the excavation process can be described as well.
Wang et al. [6] derived the yield approach index function of Hooke-Brown criterion according to the Mohr-Coulomb yield approach index expression. Moreover, the failure degree theory was introduced, and the failure degree concept of coal was defined. Finally, the crack classification method of coal was proposed by dividing the values of the failure approach index. Yao et al. [7] comprehensively and intuitively assessed the stability degree of rock and soil mass in various regions around the shield tunnel running through Hefei Metro Line 1 using the failure approach index theory and quantitatively defined the range of damage area and disturbance area. Zhou et al. [8] evaluated the similarities between the scale model of slope and the prototype failure through the failure approach index theory and proved that the inclined loading method can properly reproduce the occurrence process of landslide. On account of the improved failure approach index, Zhang et al. [9] made a quantitative analysis and assessed the influence degree of the water storage and rainfall in the Three Gorges Reservoir and their combined effects on the bank slope stability of Qianjiangping. Liu et al. [10] employed the failure approach index to quantitatively describe the risk and failure state of surrounding rock and soil mass during the construction of pipe jacking.
Evidently, the application of the failure approach index provides considerable convenience for the stability analysis of geotechnical engineering. However, the formulae for calculating the plastic shear strain and the ultimate plastic shear strain among the failure approach index calculations are complicated and difficult to solve. Meanwhile, the failure approach index represents the stability evaluation theory of surrounding rock based on the Mohr-Coulomb strain softening model, and the simple two-dimensional plastic shear strain γ p is employed in this calculation formula. With regard to the finite difference software FLAC 3D , the Mohr-Coulomb strain softening model adopts the three-dimensional plastic shear strain increment Δε ps to calculate the plastic shear strain, which brings crucial difficulties to the three-dimensional numerical simulation. Thereafter, on the basis of prior studies [11,12], the Mohr-Coulomb strength theory denotes the mathematical expression under the conditions of the criterion of positive compression and negative tension and the sequence of the principal stress σ 1 � σ 2 � σ 3 . Nevertheless, with respect to the 3D numerical simulation software, the mechanical symbol convention of positive tension and negative compression and the sequence of the principal stress σ 1 � σ 2 � σ 3 are adopted. For instance, the corresponding function of the Mohr-Coulomb yield approach index and the expression of the ultimate plastic shear strain in the calculation formula of the failure approach index cannot be found in FLAC 3D [13] and the discrete element software PFC 3D [14].
In this study, the expression of Mohr-Coulomb yield approach index is derived under the conditions of the criterion of positive tension and negative compression and the sequence of the principal stress σ 1 � σ 2 � σ 3 . Under the condition of the principal strain sequence ε 1 � ε 2 � ε 3 , the analytical expression of the ultimate elastic shear strain and the calculation formula of the ultimate elastic-plastic shear strain are derived so that the simple two-dimensional expression of the ultimate plastic shear strain is obtained. Subsequently, the relationship between the three-dimensional ultimate plastic shear strain and the two-dimensional ultimate plastic shear strain in the Mohr-Coulomb strain softening model are established to determine the critical parameters for the three-dimensional numerical simulation. Eventually, the correctness and reliability of the ultimate strain analysis method are verified by numerical calculation comparison and engineering examples. Therefore, the research findings have developed the theory of the failure approach index, which provides the analytical method and technical reference for the stability analysis of surrounding rock in underground engineering.

Mohr-Coulomb yield approach index
Sign convention: the positive and negative values of the normal stress refer to the tensile and compressive stresses, respectively (i.e. positive tension and negative compression), and the principal stress order is σ 1 � σ 2 � σ 3 . Hence, under the principal stress order of σ 1 � σ 2 � σ 3 , the Mohr-Coulomb yield criterion determined by the principal stress can be expressed as: where σ 1 and σ 3 represent the minimum and maximum principal stress, respectively; c signifies the cohesion, and φ denotes the internal friction angle. Under the conditions of the given principal stress order, the position of the point P in the plane π from the principal stress space is illustrated in Fig 1. According to the results presented in Fig 2, the projections of stress on the plane π on the x-axis and y-axis can be determined as: Upon establishing the polar coordinate system on the plane π, the expression of the stress Lord angle is obtained as: where θ σ denotes the stress Lord angle, namely, the included angle between the stress PQ on the plane π and the vertical line of the principal axis O 0 s 0 2 . Since the magnitude of the principal stress is irrelevant to the choice of the coordinate system, the expressions of the first invariant of the stress tensor I 1 and the second invariant of the stress deviator J 2 can be represented as: Combining Eqs (6), (7) and (8), the expression of the principal stress function, which is signified by the mean principal stress p, the generalized shear stress q and the stress Lord angle, is obtained as: By substituting the Eq (9) into the Eq (3), which is the yield criterion, denoted by the mean principal stress, the generalized shear stress, and the stress Lord angle, is calculated as: The spatial position relation between the stress point and the yield surface in the principal stress space is depicted in Fig 3. Moreover, the geometric relationship between the stress point

PLOS ONE
Failure approach index theory of strain limit analysis P and most stable reference point Q on the meridional plane is illustrated in the left side of Fig  3. Thereinto, the horizontal axis represents the isocline line, and the shear stress on the plane π can be expressed using the vertical coordinate τ π as: In the principal stress space, the stress state of a point P can be denoted as the point in the meridional plane and plane π. As demonstrated in Fig 3, the line segment CG stands for the meridian on the yield plane, and the coordinate values of the points Q, P, and C refer to (σ π , 0), (σ π , τ π ), and (s p ; t 0 p ), respectively. Since the point C satisfies the Eq (10) concerning the Mohr-Coulomb yield criterion, which is located on the yield surface. Accordingly, the following expression is obtained: With respect to the triangle relationship implied in Fig 3, ΔPAC and ΔQQ 0 C are similar to each other. Hence, the yield approach index can be denoted as: By substituting the Eqs (11) and (12) into the Eq (13), the expression of the yield approach index under the Mohr-Coulomb yield criterion is attained as: Where: YAI refers to the yield approach index.

Ultimate strain analysis
Rather than considering the deformation and displacement of materials, by solely investigating the material failure of rock and soil masses, the calculation process can be simplified by sidestepping the complex elastic-plastic constitutive models, and the ultimate analysis method of ideal elastic-plastic models can be adopted [15]. In the stress-strain curve of an ideal elasticplastic material, if it is expressed by stress, the stress is constant and the strain increases infinitely. In this case, the yield and failure criteria are identical, so it is difficult to distinguish the yield and failure state. In terms of the strain expression, the initial yield of the material is defined upon attaining the yield state, which has the ultimate elastic shear strain γ y . At this moment, the material will not be damaged. Nevertheless, with the development of plasticity, the material will enter a state of failure. Thus, the ultimate elastic-plastic shear strain γ f is attained, i.e., the material encounters the state of the failure strain [16].

Ultimate elastic shear strain
According to the generalized Hooke's law, where E denotes the elastic modulus, μ represents the Poisson's ratio, and ε m = (ε 1 + ε 2 + ε 3 )/ 3, which symbolizes the average normal strain. By substituting the Eq (15) into the Eq (3), the Mohr-Coulomb yield criterion expressed by the principal strain under the principal strain condition of ε 1 � ε 2 � ε 3 can be obtained as: When materials comprising rock and soil masses attain the ultimate elastic state under triaxial stress, their ultimate elastic principal strain should satisfy the following equation: where ε 1y and ε 3y denote the minimum and maximum ultimate elastic principal strain, respectively. By substituting the Eq (17) into the Eq (16), the expression of the ultimate elastic principal strain can be realized as: On the basis of the condition of the unidirectional stress, the Eq (18) can be simplified as follows: With respect to the derivation of the theoretical model, considering the simplicity of the process and calculation, the ultimate elastic shear strain is expressed as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi

Ultimate elastic-plastic shear strain
The strain space, a three-dimensional space comprising three main strains, is used to denote the strain states. As depicted in Fig 4, under the principal strain condition of ε 1 � ε 2 � ε 3 , a fixed point in the strain space corresponds to a certain strain state. By virtue of considering the rectangular coordinate O xy on the strain plane π and ensuring that the y-axis direction is consistent with the projection of the axis ε 2 on the plane π, as illustrated in Fig 5, the following expression is derived: Considering the polar coordinate of r ε , θ ε on the strain plane π, the following representative expression are derived as: where r ε depicts the strain vector diameter, and θ ε indicates the strain Lord angle. Under the condition of principal strain order of ε 1 � ε 2 � ε 3 , the principal strain function can be expressed via the variables r ε , θ ε , and ε m as: On the deviator strain plane, the relationship among the principal strain, the shear strain, and the strain Lord angle is expressed as: Under the condition of principal strain of ε 1 � ε 2 � ε 3 , given the average normal strain ε m = (ε 1 + ε 2 + ε 3 )/3, substituting this expression into Eq (25) yields the expression of the maximum shear strain γ max as: In terms of the unidirectional stress, the isotropic material, and the constant Poisson's ratio, the lateral strain of the material will be identical. In this case, the strain Lord angle θ ε = 30˚, which is substituted into the Eq (26), and the relationship between the ultimate elastic-plastic shear strain of the material and the ultimate elastic-plastic principal strain can be obtained as: where γ f represents the ultimate elastic-plastic shear strain, and ε 3f and ε 1f symbolize the ultimate strain values of the elastic-plastic principal strain of ε 3 and ε 1 . The numerical solution of the ultimate elastic-plastic shear strain can be obtained via the numerical method of the ultimate strain analysis. Upon implementing the numerical ultimate analysis method, a point failure will occur when a certain element in the material model attains the state of the ultimate plastic shear strain limit. Furthermore, if the shear strain of the material model at all points on the penetrating failure plane is greater than the ultimate plastic shear strain, the overall failure state of the material model will emerge, and the penetration of the ultimate strain zone represents the sufficient and necessary condition for the overall failure state of the material model [17].

Ultimate plastic shear strain
The failure and yield surface of the Mohr-Coulomb possess the same shape and size, an unequal-angle hexagonal cone, with an unequal-angle hexagon on the partial plane, and the distance between the center of failure surface and center of yield surface is equal to g p f [18]. According to the ideal elastic-plastic model, the expression of the ultimate plastic shear strain can be represented as: In relation to the Mohr-Coulomb's strain-softening model, the plastic parameters are introduced to account for the dependence of the loss of rock strength on plastic deformation. Generally, common plastic parameters are in two forms: the form of the internal variable, namely, γ p , and the the form of the plastic strain increment, i.e. Δε ps .
Under the condition of principal strain order of ε 1 � ε 2 � ε 3 , the two-dimensional plastic shear strain can be expressed as [19]: where γ p represents the two-dimensional plastic shear strain, and ε ps 3 and ε ps 1 signify the maximum and the minimum plastic principal strain, respectively.
In the Mohr-Coulomb's strain-softening model, the three-dimensional plastic shear strain is adopted [19,20], whose incremental form is presented as follows: Dε ps ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Dε ps where Dε ps m ¼ Dε ps 1 þ Dε ps 3 ð

Þ=3: ð31Þ
As for the finite difference software, FLAC 3D , the strain-softening model based on the Mohr-Coulomb criterion adopts the three-dimensional plastic shear strain increment. Moreover, establishing a relationship between the internal variable, γ p , and the three-dimensional plastic shear strain, Δε ps , in the strain-softening model is essential for embedding the failue approach index formulation into the FLAC 3D software for numerical calculation. The Mohr-Coulomb plastic potential function and yield function have similar expressions. Under the condition of principal stress of σ 1 � σ 2 � σ 3 , the shear plastic potential function can be expressed as where g s represents the shear plastic potential function, and ψ denotes the dilatancy angle.
In the Mohr-Coulomb's strain-softening model, the shear plastic potential function corresponds to the law of uncorrelated flow, which is expressed as: where λ s denotes the plastic factor.
According to the Eq (34), the relationship between the principal strain increments Dε ps 3 and Dε ps 1 is established, which can be signified as: In accordance with the Eqs (31) and (35), the Eq (30) can be simplified to the following: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi In case of a constant dilatancy angle, the relationship between the three-dimensional plastic shear strain, ε ps , and the two-dimensional plastic shear strain, γ p , in Mohr-Coulomb's strainsoftening model can be denoted as: Based on the Eq (35), the Eq (38) can be rewritten as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi With regard to the relationship between the Eqs (38) and (39), upon combining with the Eq (29), the three-dimensional plastic shear strain can be represented as: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi Assuming that the dilatancy angle ψ = 0, the Eq (40) can be simplified as: In line with the ultimate strain state, the three-dimensional ultimate plastic shear strain in the Mohr-Coulomb' strain-softening model is expressed as:

Solution and verification of ultimate strain
In this study, the numerical ultimate analysis of C20 cube concrete specimen is carried out by FLAC 3D . The geometric size of the calculation model is 150 mm×150 mm×150 mm. Each side is divided into 19 grids, with a total of 6859 elements and 8000 nodes. Three direction constraints are applied to the bottom surface of the model, and uniform vertical load downward is applied to the top surface without considering the friction of the top surface. According to the failure form of the cube concrete specimen, there are 12 key recording points (units) in the calculation model, as shown in Fig 6. In addition, the physical and mechanical parameters of C20 concrete are listed in Table 1.
Through the numerical simulation analysis, several criterion methods for failure have been proposed [22,23], and the method of iterative non-convergence for failure estimation has been employed in this study. In the calculation process, the axial load of the model is continuously increased until the model approaches the state of ultimate failure, at which point the load signifies the ultimate load. Fig 7(a) depicts the maximum unbalanced force curve for a load of 20.03 MPa. Evidently, the latter section of the convergence curve presents a horizontal straight line and tends to be zero, indicating that the calculation converges and the model will not be damaged. Similarly, the maximum unbalanced force curve for a load of 20.04 MPa is depicted in Fig 7(b): the latter part of the curve fluctuates continuosuly. In particular, observing the convergance trend during the calculation is a complicated task, implying that the model will be damaged.
Therefore, the ultimate load onto the concrete cube specimen is 20.03 MPa, which is consistent with the strength grade of the C20-grade concrete, and this result indicates that the given shear strength parameters of concrete are accurate and reliable.
Generally, the elastic strain of concrete is linearly related to the stress level and can be obtained directly using Hooke's law. Under the condition of principal strain sequence of ε 1 � ε 2 � ε 3 , the ultimate elastic principal strain and ultimate elastic shear strain of concrete can be calculated using the Eqs (19) and (20), respectively, and the values of the calculated ultimate strain are enumerated in Table 2.
Concrete is regarded as an ideal elastoplastic material, and the ultimate strain values of 12 elements in the model are calculated under the ultimate load. During the numerical calculation, the development of the strain of unit 5, which is onto the middle edge of the top surface of the model, was evident. When the model is cycled to the default maximum unbalanced force ratio, the strain of the unit 5 approaches the maximum level and failure will inevitably occur. The curve of the axial load and principal strain are exhibited in Fig 8. As the unit is located at the central portion of the top surface of the cube, the failure of the unit will inevitably lead to overall failure of the cement-based concrete cube specimen. Notably, this element strain represents the ultimate strain of the C20-grade concrete. Under the condition of principal strain order of ε 1 � ε 2 � ε 3 , the ultimate elastoplastic strains of concrete are listed in Table 2.  Table 1. Physical and mechanical parameters of C20-grade concrete [21].  The ultimate plastic shear strain of concrete is calculated according to the Eq (28), whose value is displayed in Table 2.

Concrete Density ρ/(kg/m 3 ) Elastic Modulus E/(GPa) Possion's ratio μ Cohesion c/(MPa) Angle of internal friction φ/(˚) Tensile strength σ t /(MPa
The ultimate strain values of concrete, as summarized in the monograph "Concrete Structure Design Principle", are listed in Table 3 [24]. On the basis of the condition of principal strain order of ε 1 � ε 2 � ε 3 , the ultimate strain values of concrete calculated by the aforementioned formulas are consistent with those in Table 3, which verifies the correctness of the ultimate strain analysis and reliability of the numerical ultimate analysis method.

Application of failure approach index theory of ultimate strain analysis to stability analysis of coal pillar
The 9-10 coal seam of Linfen Tianyu Hengjin Coal Industry Co., Ltd., Anhui North Coal and Electricity Group, is commercially mined; the average thickness of the mine is 5.6 m, and the dip angle is 3-6º, which corresponds to that of a near-horizontal type. Moreover, the ground elevation ranges from +1397 m to +1424 m, and the mining upper-and lower-limit elevation is +1043.32 m and +983.34m, respectively. During the extraction process, the width of the coal pillar in the remaining section is approximately 30 m as yet, which is the traditional empirical width. Considering the geological conditions of the shallow buried depth, thick coal seam, and hard roof, this wide coal pillar will cause the waste of resources, and thus, the stability of narrow coal pillar should be investigated from a roadway protection perspective to reduce the wastage of coal resources. Hence, the determination of a reasonable width for coal pillars is of great significance to the coal mining industry in this area.
In this study, the lateral driving roadway of 9101 goaf in Xiyi mining area is considered as the research object to investigate the stability of the section coal pillar. The length of the 9101 working face is 120 m, and the main roof is constituted of K2 limestone with an average thickness of 13.09 m. Additionally, the immediate floor is mudstone, and the local face comprises sandy mudstone with an average thickness of 1.82 m. Moreover, the main floor is mainly comprises aluminum mudstone and limestone with an average thickness of 12.71 m. The tunnel section (length: 4 m; height: 2.5 m) represents a rectangular form. The physical and mechanical parameters of each rock layer are summarized in Table 4.
By means of the Eq (20), the ultimate elastic shear strain is calculated to be 0.16%. Obtained via the numerical analysis of ultimate strain, the 12 key units of the elastic and plastic shear strain curve in the calculation model are exhibited in Fig 9. Notably, the elastic and plastic shear strain of the unit 4 is developing promptly, which is located in the upper specimen. In this case, the failure of the unit is bound to result in the overall failure of the specimen. Thus, the strain of this unit denotes the ultimate elastoplastic shear strain, whose value is equal to 1.3%. Accordingly, the ultimate plastic shear strain of the coal pillar is 1.15%. According to the Eq (42), the three-dimensional plastic shear strain value of the Mohr-Coulomb's strain-softening model is 0.0057, which is adopted in the FLAC 3D software for the numerical simulation. Ultimately, the elastoplastic quantitative value of the section coal pillar can be obtained using the failure approach index theory, thereby yielding an updated method for determining the reasonable width of the section coal pillar. On account of the aforementioned ultimate strain analysis, the ultimate plastic shear strain value of the coal pillar is determined, which provides critical values for the calculation of failure approach index. The numerical simulation flow chart of the failure approach index theory of the ultimate strain analysis in section coal pillar stability analysis is displayed in Fig 10. Evidently, there are 119611, 118586, and 125592 nodes, and 693631, 688130, and 729879 units of the coal pillar model for widths of 5, 7, and 9 m, respectively. In addition, the boundary conditions of the coal pillar model signify that the vertical downward uniform load is imposed on the top surface, and the other five surfaces are fixed.
In the coal pillar bearing core area, the failure approach index of the three schemes is 1.17, 1.1, and 0.9, as demonstrated in Fig 11. According to the elastic kernel theory of the coal pillars [26] and the failure approach index partition [5], the width of the coal pillar for the failure approach index of 0.9 is considered the reasonable width of the coal pillar in the reserved section. Therefore, the width of the coal pillars in Xiyi mining area is approximately 9 m.
To verify the own stability of width with 9 m of coal pillar in Xiyi mining area, a monitoring station on Section I of surface displacement and anchor axial force is arranged at a distance of 5 m from the roadway head during the excavation of the mining roadway. By virtue of monitoring the force change of coal pillar and observing the deformation of the roof and floor of the roadway and the two sides of the roadway, the stability of the coal pillar is fed back, and the rationality of the width of the coal pillar is judged. As indicated by Fig 12, the axial force of the bolt increases rapidly when the roadway is pushed about 20 m, and the change rate of the axial force of the bolt slows down significantly when the roadway is pushed about 100 m from the head. At 150 m from the head, the axial force of the bolt gradually becomes stable. As can be observed from Fig 13, the cumulative displacement of the top and bottom plate of the measuring station I is 155 mm, that of the solid coal side is 60 mm, and that of the coal pillar side is 75 mm, which is consistent with the numerical displacement nebulogram of the coal pillar side. With the advance of roadway driving, the influence of roadway driving on the measuring station I gradually decreases, whose deformation and deformation rate tend to be stable, the  roadway rock pressure is not proved to be apparent, and the surrounding rock stability of the roadway performs preferable. Consequently, the research results of axial force and surface displacement of anchor bolt in station I indicate that it is reasonable to reserve the section with coal pillar of width of 9 m.

Discussion
Currently, the rapid development of the ultimate analysis and numerical ultimate analysis methods have provided a convenient way to solve the ultimate strain of geotechnical materials, whose reliability and adaptability have also been demonstrated [27,28]. The failure approach index theory of the ultimate strain analysis proposed in this paper offers a simple expression of the ultimate plastic shear strain and solves the problem of the value of the ultimate plastic shear strain through the three-dimensional numerical simulation of the failure approach index, which provides convenient conditions for the broader application of the failure approach index and improves and develops the theory of the failure approach index.
Herein, the ultimate shear strain refers to the expression obtained under the condition of unidirectional force. If considering the triaxial stress or anisotropy of rock and soil materials, the methodology in this paper can be improved. However, the expression will evidently be more complicated, and the solution will become more difficult. Meanwhile, the failure approach index theory of the ultimate strain analysis put forward in this study can be applied to rock underground engineering, whose application scope needs to be further expanded. It can also be applied to pipeline engineering, subway tunnel engineering, deep underground disposal storage of nuclear waste and other construction projects for extending the application scope and examining its applicability.

Conclusions
1. Notably, the Mohr-Coulomb yield criterion in the principal stress space is derived under the conditions of positive tension and negative compression and the sequence of the principal stress σ 1 � σ 2 � σ 3 . Moreover, the expression of the Mohr-Coulomb yield approach index is obtained according to the spatial geometric relation.
2. Under the principal strain sequence ε 1 � ε 2 � ε 3 and the condition of the unidirectional force, the analytical expression of the ultimate elastic shear strain is derived, and the relationship between the ultimate elastic-plastic shear strain and the principal strain is acquired. Then, the two-dimensional expression of the ultimate plastic shear strain is obtained. Furthermore, the relationship between the three-dimensional plastic shear strain and the two-dimensional plastic shear strain in the Mohr-Coulomb strain softening model is established so that the calculation parameters of the three-dimensional ultimate plastic shear strain in the Mohr-Coulomb strain softening model can be determined and the failure approach index theory can be improved.
3. The ultimate strain analysis of specimen of concrete C20 is carried out through the derived theoretical formulae, and the obtained ultimate strain values are found to be consistent with the previous research findings, which verifies the correctness and reliability of the ultimate strain analysis method.
4. The failure approach index theory of the ultimate strain analysis is applied to the quantitative elastic-plastic analysis of the section coal pillar in the section of Hengjin coal industry. The reasonable retainment width of the section coal pillar is calculated as 9 m, which is shown to be rational by using the feedback of the failure approach index theory from the monitoring results of roadway mine pressure.
Supporting information S1 Table. Table data